package optimizers.tools;

//author: nicolas.bredeche(at)isir.upmc.fr
//created 2013-2-8

public class GaussianRandom {  // box-muller implementation comes from the internet

	// internal buffer (Box-Muller generates random numbers two by two)
	boolean 	_buffered = false;
	double 		_bufferValue;
	
	// random number will be generated according to N(mu,sigma^2)
	double 		_mu;
	double		_sigma;

	public GaussianRandom ( )
	{
		_mu = 0;
		_sigma = 1;
	}

	public GaussianRandom ( double mu, double sigma )
	{
		_mu = mu;
		_sigma = sigma;
	}
	
	// basic form of Box-Muller transformation (constant time but stability issue when x1 close to zero)
	private double[] generateNumbersSlow( double x1 , double x2 )
	{
		double res[] = new double [2];
	    res[0] = Math.sqrt( - 2 * Math.log(x1) ) * Math.cos( 2 * Math.PI * x2 );
	    res[1] = Math.sqrt( - 2 * Math.log(x1) ) * Math.sin( 2 * Math.PI * x2 );
	    return res;
	}

	// polar form of Box-Muller transformation (faster, reliable)
	private double[] generateNumbersFast( double x1 , double x2 )
	{
	    double w;
		double res[] = new double [2];
	
	    do {
	            x1 = 2.0 * Math.random() - 1.0;
	            x2 = 2.0 * Math.random() - 1.0;
	            w = x1 * x1 + x2 * x2;
	    } while ( w >= 1.0 );
	
	    w = Math.sqrt( (-2.0 * Math.log( w ) ) / w );
	    
		// res[0] = mu - sigmaT1 et res[1] = mu - sigmaT2 suivent la loi normale N(mu,sigma^2)
	    res[0] = this._mu - this._sigma * ( x1 * w );
	    res[1] = this._mu - this._sigma * ( x2 * w );
	    
	    // System.out.println( "x1= "+ res[0] + " ; x2= "+ res[1] );
	    
	    return res;
	}

	/**
	 * Gaussian number generator
	 * @return random value between 0 and 1 wrt. to gaussian law
	 */
	public double random()
	{
		if ( this._buffered == true )
		{
			this._buffered = false;
			return this._bufferValue; 
		}
		else
		{
			double res[] = generateNumbersFast(Math.random(),Math.random());
			this._bufferValue = res[1];
			this._buffered = true;
			return res[0];
		}
	}
	
	// accessing methods 
	
	public void setParameters(double __mu, double __sigma) 
	{ 
		this._mu = __mu; 
		this._sigma = __sigma; 
		this._buffered = false; 
	}
	public double getMu() { return this._mu; }
	public double getSigma() { return this._sigma; }
	

	// debug purpose
	
	public static void main(String[] args) {
        
		double startTime = System.currentTimeMillis();
		
		GaussianRandom gaussianrandom = new GaussianRandom(0,0.33);
		GaussianRandom gaussianrandom2 = new GaussianRandom(0.5,0.5);
		
        for ( int i = 0 ; i != 100 ; i++ )
        	System.out.println("1,"+gaussianrandom.random()+",2,"+gaussianrandom2.random());
        
        //Display.info("\nTerminated ("+ ((System.currentTimeMillis()-startTime)/1000) +"s elapsed).");
    
	}
}


/*

# create a graph from an "exp.data" file
# example for plotting data from gaussian law random number generator
#
# syntaxe:  gnuplot plotGauss.gp
#
# no_iteration learning_error generalization_error
# an output file "exp.eps" is created

set xlabel 'x'
set ylabel 'y'
set title 'random number with gaussian law'
set key left bottom
set key box
set datafile separator ","

set yrange[-3:3]
set xrange[0.8:2.2]

plot 'expGauss.dat' using 1:($2) title "N(0,1)" with points, 'expGauss.dat' using 3:($4) title "N(0.5,0.5)" with points

# Decommenter ce qui suit pour generer un fichier EPS en sortie

#set term post eps "Times-Roman" 8
#set size 5./10., 3./7.
#set output 'exp.eps'
#replot
#set term X11

# OU ALORS:
set term postscript eps enhanced monochrome
set output 'exp.eps'
replot

pause -1

*/